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Abstract 

The evolution of spherical single-mass star clusters driven by two-body relaxation was 
followed beyond core collapse by numerically solving the orbit-averaged Fokker-Planck equa- 
tion in energy-angular momentum space. The heating effect by three-body binaries was 
incorporated in the Fokker-Planck models. The development of velocity anisotropy after 
the core collapse is discussed in detail. The anisotropy in the outer regions continues to 
increase slowly after the collapse. In clusters comprising a relatively small number of stars 
(N < 10000), the post-collapse expansion is nearly self-similar and the anisotropy at each 
of inner Lagrangian radii is nearly constant. In clusters comprising a larger number of stars 
(N > 10000), gravothermal oscillations occur and the anisotropy at each of the inner radii 
oscillates with the core oscillations. There is no qualitative difference in the nature of the 
gravothermal oscillations between the isotropic and anisotropic Fokker-Planck models. 

Key words: Clusters: globular — Fokker-Planck equation — Numerical methods - 
Stars: stellar dynamics 



1 Introduction 



The dynamical evolution of globular clusters has been extensively investigated [see Spitzer 
(1987) for a review]. In many of those investigations, Cohn's (1980) direct integration scheme 
for the Fokker-Planck (hereafter FP) equation was used as a main numerical tool. Recently, 
many studies have been made particularly to reveal the realistic evolution of globular clusters 
and to compare the theoretical models with observations. In such studies various effects were 
incorporated in numerical simulations: the stellar-mass spectrum, binaries, the galactic tidal 
field, stellar evolution, etc. (e.g., Chernoff, Weinberg 1990; Drukier 1995). On the other hand, 
the anisotropy of the velocity distribution was almost always neglected, although it is obvi- 
ous that the anisotropy develops at least in the halo (it is expected that the radial velocity 
dispersion exceeds the tangential one). This neglect was mainly due to a numerical difficulty 
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involving anisotropic FP models. The evolution of anisotropic clusters can be described by a 
two-dimensional (2D) orbit-averaged FP equation in energy-angular momentum space (Cohn 
1979), while the evolution of isotropic ones can be described by a one-dimensional (ID) orbit- 
averaged FP equation in energy space (Cohn 1980). The direct integration code for the 2D FP 
equation had a numerical problem in that the energy conservation was insufficient to continue 
the run beyond a factor of 10 3 increase in the central density (Cohn 1979; Cohn 1985). On the 
other hand, an integration of the ID FP equation can be performed with much higher numer- 
ical accuracy (Cohn 1980), largely due to the adoption of the Chang-Cooper finite differencing 
scheme (Chang, Cooper 1970). 

Recently, Takahashi (1995, hereafter Paper I) has developed a numerical method for solving 
the 2D FP equation. The method is essentially the same as Cohn's (1979) method. A main 
difference between the two methods exists concerning discretization schemes of the FP equation. 
Cohn (1979) used a finite-difference scheme in which simple centered-differencing was adopted 
for the spatial discretization. Cohn (1985) reported that he investigated several heuristic gener- 
alizations of the Chang-Cooper scheme, and that all of these improved the energy conservation, 
though the details of these schemes were not explained. In Paper I, two different discretization 
schemes were employed: one was a finite-difference scheme where the Chang-Cooper scheme is 
simply applied for only the energy direction; the other was the finite-element scheme, where 
the test and weight functions implied by the generalized variational principle (Inagaki, Lynden- 
Bell 1990; Takahashi 1993) are used. Using those schemes, the gravothermal core collapse was 
followed until the central density increased by a factor of 10 14 with a 1% numerical accuracy 
concerning the total-energy conservation. This was a big advance compared with previous cal- 
culations; the central density growth factors in the calculations of Cohn (1979) and Cohn (1985) 
were 10 3 and 10 6 , respectively. It should be noted that a numerical error originates not only in 
the integration of the FP equation, but also in other calculation procedures, e.g., the calculation 
of the diffusion coefficients and the potential-recalculation steps. It should also be noted that 
2D FP calculations require a rather large computational time (see section |3|). Thus, ten years 
ago it was not easy to perform such calculations as those which we present here. Besides the 
FP models, anisotropic gaseous models of star clusters have recently been successfully applied 
(e.g., Giersz, Spurzem 1994; Spurzem 1996). 

In Paper I, the pre-collapse evolution of single-mass clusters was studied. In particular, 
Paper I revealed the evolution during self-similar phases of core collapse in anisotropic clusters. 
The density profile left outside the collapsing core is the same as that in isotropic clusters; i.e. 
p oc r~ 2 23 . In the self-similar regions, a slight velocity anisotropy exists: i.e. of /of = 0.92, 
where a r and at are the one-dimensional radial and tangential velocity dispersions, respectively. 
The core collapse rate, £ = t T (0)dln p(0) / dt, where i r (0) and p(0) are the central relaxation time 
and density, is £ = 2.9 x 10~ 3 , which is 19% smaller than the value of £ = 3.6 x 10~ 3 for an 
isotropic model. That is, the core collapse proceeds slightly more slowly in the anisotropic model 
than in the isotropic model. When the initial model is Plummer's model, the collapse occurs 
at time 17.6i r h,i in the anisotropic model and 15.6i r h,i in the isotropic model, where i r h,i is the 
initial half-mass relaxation time. The halo soon becomes dominated by radial orbits, even if the 
velocity distribution is initially isotropic everywhere. The ratio of the radial velocity dispersion 
to the tangential one increases monotonically as the radius increases. 

Following Paper I, this paper examines the post-collapse evolution of single-mass clusters. 
The effect of three-body binaries is incorporated into FP models as a heat source. We are 
particularly interested in the development of the anisotropy in the halo after the core collapse. 
Does the anisotropy continue to increase even after the collapse, or come to be constant ? We 
are also interested in whether there are any differences concerning the nature of gravothermal 
oscillations between isotropic and anisotropic models. In section 2, the models and numerical 
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methods are described. In section 3, the calculation details are described and the numerical 
accuracy is discussed. Section 4 presents the results of the calculations. The conclusions and 
discussion are given in section 5. 

2 The Models and Methods 

2. 1 Fokker-Planck Models 

We consider the collisional evolution of spherical single-mass star clusters in dynamical equi- 
librium. In such clusters the distribution function of stars, /, is a function of the energy per unit 
mass E, the modulus of the angular momentum per unit mass J, and time t; i.e. / = f(E, J, t). 
The evolution of / due to the two-body relaxation can be described by an orbit-averaged FP 
equation in (E, J)-space (Cohn 1979). Numerical integration of the FP equation is performed 
in the same manner as in Paper I, but a binary heating term is included in the equation. 

For problems concerning post-collapse evolution, we must specify the number of stars in the 
cluster, N, and the numerical constant, fj,, in the Coulomb logarithm ln(/xiV) (e.g. Spitzer 1987, 
p30). In all of the calculations we adopted /x = 0.11, which was obtained by Giersz and Heggie 
(1994a) for the pre-collapse evolution of single-mass clusters. We note that Giersz and Heggie 
(1994b) found a smaller value of \i for the post-collapse evolution (their best value was [i = 0.035). 
However, we fixed the value of \i throughout all evolutionary phases. A small difference in fi 
does not seriously affect the nature of cluster evolution. Although the determination of an 
appropriate value of p is an interesting subject in collisional stellar dynamics, it was beyond the 
scope of this study. A future careful comparison between iV-body, gaseous, and FP models may 
give further information concerning the Coulomb logarithm. 

2. 2 Three-Body Binary Heating 

The three-body binary heating rate per unit mass is given by 

E h = C h G 5 m 3 p 2 a- 7 (1) 

(Hut 1985), where m is the stellar mass, p the mass density, a the one-dimensional velocity 
dispersion, and Cb a numerical coefficient. In this paper we choose the standard value of 
Cb = 90. The local heating rate ([[]) is orbit- averaged as 

f r& dr ■ I f r& dr 

(-E'b)orb = / E h / , (2) 

Jr p V v I J rp V v 

where v r = {2[<^>(r) — E] — J 2 /r 2 } 1 ^ 2 is the radial velocity of a star of energy E and angular 
momentum J at radius r, and r p and r a are the pericenter and apocenter radii of the star, 
respectively. The orbit-averaged heating rate (g) is added to the usual first-order diffusion 
coefficient (AE) OT ^, (cf. Cohn et al. 1989). Furthermore, we assume that the scattering by 
binaries does not produce the net changes of the scaled angular momentum R, i.e. (i?b)orb = 0. 
Here, R is defined as R = J 2 /J^(E), where J C (E) is the angular momentum of a circular orbit 
of energy E. 

3 Numerical Calculations 

Plummer's model (e.g. Spitzer 1987, pl3) was chosen as the initial cluster model, where 
the velocity distribution is isotropic everywhere. Test calculations were carried out using both 
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the finite-difference and finite-element codes described in Paper I. In calculations of the pre- 
collapse evolution, the two codes achieved similar numerical accuracy concerning the total energy 
and mass conservation, and the results obtained using them were generally in good agreement 
(Paper I). In the present calculations of the post-collapse evolution, the results obtained by the 
two codes were also generally in good agreement. However, the numerical error in the energy 
conservation was considerably larger in the case of the finite-element code. 

We note that the total energy of the cluster cannot be conserved in these calculations, but 
should increase, because the binary heating is taken into account. The total energy should 
increase by the amount of energy input. To check the numerical error, the energy input was 
recorded during the calculations. The energy input at each time step may be calculated by 
integrating the product of the binary heating rate [equation @] and the distribution function 
over energy-angular momentum space. The cumulative energy input is summed up over the 
time steps of the run. There must be some degree of inaccuracy inherent in this estimation for 
the energy input. However, this estimated energy input and the actual energy increase resulting 
from the FP-integration should be in agreement within some degree of numerical accuracy; the 
degree of the agreement becomes better as the mesh sizes and the time step size become smaller. 
When we estimated the numerical error in the energy conservation, we assumed that the energy 
input estimated as above was exact. 

For example, at the end of the calculation for N = 5000 with the finite-difference code 
(see figure la), the relative energy error, which is defined as the ratio of the amount of the 
energy change due to numerical error to the initial total energy, was about 2%; at the end of a 
corresponding calculation with the finite-element code, the relative energy error was about 12%. 
There was a systematic energy drift during calculations of the post-collapse evolution in both 
the finite-difference and finite-element codes. This error arose mainly from the FP-calculation 
steps. In fact, in one FP step, the actual energy increase was always slightly smaller than the 
expected increase due to binary heating. The degree of this disagreement was much larger in 
the finite-element code than in the finite-difference code. An energy error arose also from the 
Poisson-calculation steps. However, the sign of the energy change in one Poisson step was nearly 
random, and the sum of the changes was small. Therefore, the energy error stemming from the 
Poisson steps does not contribute very much to a cumulative error. 

The reason why the accuracy of the finite-element code for the energy conservation is not 
very good for the post-collapse calculations is not yet clear. One way to improve the accuracy 
is to increase the mesh numbers, especially for the energy. In fact, the energy error decreased 
as the energy- mesh number increased, although the accuracy of the finite-difference code was 
better with the same mesh number. Another promising way to improve the accuracy of the 
finite-element scheme is to use higher-order basis functions (see appendix 2 of Paper I). In the 
present code two-dimensional piecewise bilinear polynomials are used as the basis functions. 
The use of higher-order basis functions, however, introduces rather complicated computational 
procedures, and, as a result, a larger computational time. In addition, we can obtain reasonably 
good accuracy with the finite-difference code. Thus, we did not try higher-order basis functions 
in the present work. 

As a result, we adopted the finite-difference code for the calculations presented in this paper 
because of its higher numerical accuracy in energy conservation. The results of calculations for 
N = 5000, 10000, 20000, and 40000 are shown in section ||. We denote the number of grid points 
in X, Y, and r by Nx, Ny, and N r , respectively. [Variables X and Y are used instead of E and 
R in the code (Paper I).] In these calculations, we set Nx = 151, Ny = 35, and N r = 91. The 
radial grid was constructed between 10 -7 r$ and 10 2 Vq, where tq is a length scale parameter of 
Plummer's model. We carried out several test calculations with other sets of grid numbers, and 
confirmed that the results converged. The relative energy errors at the ends of the calculations 
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shown in figure 1 were 2.3% (N = 5000), 2.5% (N = 10000), 1.1% (N = 20000), and 0.1% 
(N = 40000). In the calculation of N = 40000, the energy error reached its maximum (0.5%) 
at the first collapse time. However, since the sign of the energy error changed with the core 
oscillations, there was some cancellation in the cumulative error. The relative mass errors were 
0.85% (N = 5000), 0.87% (N = 10000), 0.51% (N = 20000), and 0.57% (N = 40000). 

One inevitable disadvantage of the 2D FP model relative to the ID FP model is that 2D 
calculations take a much larger computational time than do ID calculations. One may expect 
that 2D calculations require about Ny times as large a computational time as do ID calculations. 
In fact, however, it can happen that the computational time of 2D calculations increases faster 
than as Ny ■ The computational time required to solve a linear matrix equation for a discretized 
FP equation is not negligible, but, rather, a few tens of a percent of the total computational 
time in 2D calculations. In ID calculations, in contrast, it is almost negligible compared with 
the total computational time, because the matrix is tridiagonal and can be inverted very easily. 
In 2D calculations, the matrix is a band matrix whose half-bandwidth is about min(iVx, Ny), 
or Ny in our cases. We can choose various direct or iterative schemes for solving the matrix 
equation. In some direct schemes for band matrices, the number of required operations varies 
as NxN Y for large Nx and Ny. A kind of conjugate gradient method (iterative method) was 
actually used in our 2D FP code. The number of operations varies as NxNy for this method. 
We found by experience that the computational time required by a 2D calculation with our 
code is about 2Ny -times larger than that required by a corresponding ID calculation (with the 
same Nx)- Most of the numerical calculations were performed on a HP 9000/715 workstation 
(at 50 MHz clock cycle). For example, 2D FP calculations for N = 5000 and 40000 (cf. figure 
1) required about 29 and 140 hours of CPU time on this machine, respectively. The total 
numbers of potential-recalculation time steps (the FP time-step size was 1/10 of the potential- 
recalculation time-step size) in these runs were 3000 and 15000, respectively, and thus the 2D 
FP code required about 34 CPU sec per step. 

4 Results 

The results are presented in standard units such that G = M = 1 and 6\ = 1/4, where G is 
the gravitational constant, M is the total mass, and E\ is the initial total binding energy. The 
time is usually measured in units of the initial half- mass relaxation time t r h,i (Spitzer 1987, p40). 

Figure 1 shows the evolution of the central density p(0) for the cases of N = 5000, 10000, 
20000, and 40000 for the ID and 2D models. For each N, the features of the evolution in the 
ID and 2D models are very similar. An apparent difference between the two models exists in 
the core collapse times. For every N, the core collapse (or bounce) occurs slightly earlier in 
the ID model than in the 2D model, as found in Paper I. Although an intuitive explanation 
for a slower collapse in the 2D (i.e. anisotropic) models is given in Paper I (see also Louis 
1990), a more convincing proof for it is desirable. There is a possibility that the slower collapse 
may be due to the numerical inaccuracy. We may be able to test this possibility simply by 
repeating the calculations with finer grids. We found that the collapse time was not affected by 
increasing the grid numbers. This fact supports the conclusion that the slower collapse rate in 
the anisotropic models is real. We also note that it is uncertain whether adopting the isotropic 
distribution function to calculate the diffusion coefficients (Paper I) has any noticeable effects 
on the collapse rate. 

The core expansion after the core collapse is stable for iV = 5000, and overstable for N = 
10000. (For N = 10000, the core expansion in the ID model is really overstable, though the 
growth of the instability is slower than in the 2D model.) For N = 20000, the core expansion is 
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unstable; the central density oscillates chaotically with a large amplitude; that is, gravothermal 
oscillations (Bettwieser, Sugimoto 1984) occur. The gravothermal oscillations also occur for 
N = 40000 with a larger amplitude than for N = 20000. Such a change in the nature of the 
post-collapse core evolution from monotonic expansion to chaotic oscillations with increasing 
./V was discussed in detail by Goodman (1987), Heggie and Ramamani (1989), Breeden et al. 
(1994), as well as Breeden and Cohn (1995), where isotropic models were used. Spurzem (1994) 
presented long-lasting gravothermal oscillations in his anisotropic gaseous model. We see no 
qualitative difference concerning the features of the gravothermal oscillations between the ID 
and 2D models. The amplitudes and periods of the oscillations, and the appearance of multiple- 
peaks in the two models are similar. This is a reasonable result, because the stability of the core 
expansion is determined by the degree of central concentration (Goodman 1987). Furthermore, 
the velocity distribution is isotropic in the core, even in anisotropic models. It is interesting 
that Spurzem (1994) suggested that gravothermal oscillations are more regular in the anisotropic 
model than in the isotropic one. 

Figure 2a shows the evolution of the Lagrangian radii containing 1, 2, 5, 10, 20, 30, 40, 50, 
75, and 90% of the total mass for N =5000. The core radius r c is also plotted; it is defined as 
n i/2 

(3) 



r c . = 



4vrG/9(0) 



(Spitzer 1987, pl6), where v m (0) is the total velocity dispersion at the center. Because of the 
difference in the collapse time, a comparison between the ID and 2D models concerning the 
evolution of the spatial structure is somewhat complicated. Thus, we plot figure 2b, where the 
time of the ID calculation is scaled so that the collapse time in the ID model should coincide with 
that in the 2D model. Concerning the evolution before the core bounce, we see small difference 
between the two models in figure 2b, except for the 90% radius. After the core bounce, the core- 
halo structure is more developed in the 2D model; that is, the 2D model has more concentrated 
inner Lagrangian radii and more extended outer radii. We note that there is no big difference 
between the two models in the evolution of the half-mass radius. The evolution of the half-mass 
radius is, roughly speaking, determined by the change in the total energy when the total mass is 
conserved. In fact, the histories of the total energy changes in the two models are similar, if the 
time is scaled as in figure 2b. In this respect, the coincidence of the evolution of the half-mass 
radius is reasonable. The effects of the development of the anisotropy on the density is apparent 
in the outer half-mass region. This is a consequence of the development of radial orbits that 
the outer Lagrangian radii are more extended in the 2D model. The more concentrated inner 
Lagrangian radii are a necessary reaction to that. However, the evolution of the core radius in 
the two models is again almost identical (if the time of the ID models is scaled). 

Figure 2b shows that the post-collapse evolution well after the core bounce seems to be self- 
similar in both the ID and 2D models; all of the Lagrangian radii as well as the core radius expand 
nearly self-similar ly. A simple argument gives a self-similar expansion law, r oc i 2 / 3 , for isolated 
clusters with no mass-loss (Henon 1965; Goodman 1984). Our 2D model as well as our ID model 
is consistent with this law. In the cases of other iV's, the evolution of the outer Lagrangian radii 
is similar to that in the case of N = 5000. When gravothermal oscillations occur, the inner 
Lagrangian radii oscillate, while the mean trend of these radii is also an expansion (cf. figure 
5). 

Figures 3a and 3b show the evolution of the anisotropy parameter A, 



A = 2- 2— | (4) 



at the 1, 2,..., and 90% Lagrangian radii, for the case of N = 5000. During the core collapse, 
A increases at every Lagrangian radius. Even in the very inner regions (e.g. at the 1 and 2% 
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radii) the anisotropy increases at advanced stages of the collapse. Just after the core bounce, 
the anisotropy at each of the inner (1-20%) Lagrangian radii decreases rapidly. This is due 
to a rapid core expansion. The core expansion is faster than the expansion of the Lagrangian 
radii located outside the core at that time. Then, the radial velocity dispersion decreases faster 
than the tangential one outside of the core, because the former is influenced more by the core 
condition. This is an exactly opposite process to that occurring during the core collapse. After 
the rapid core expansion phase, the cluster expands nearly self-similar ly (as mentioned above), 
and the anisotropy at each inner Lagrangian radius settles to roughly a constant value. 

The anisotropy at the outer Lagrangian radii continues to slowly increase after the core 
bounce. In figure 3b we can see that the curve of the anisotropy at the 90% radius flattens at 
late times. This is partly because A cannot exceed two, by definition. In any case, it is true 
that the rate of the anisotropy increase at the outer radii slows down. The development of the 
anisotropy in the outer regions is a consequence of the emergence of radial-orbit stars which have 
gained energy as the result of relaxation in the inner regions (Paper I). Therefore, we expect 
that the rate of the anisotropy increase is related to the relaxation time in the inner regions. 

Figure 4 shows the evolution of the anisotropy at the Lagrangian radii for N = 5000 as a 
function of the elapsed number of actual central relaxation times, 



(Cohn et al. 1989), where t rc is the central relaxation time. The core bounce occurs at r ~ 2400. 
This figure indicates that the anisotropy at the outer Lagrangian radii increases roughly linearly 
with r after a rapid increase at the very initial stages (r < 1000). While we see in figure 3b 
that the rate of increase of the anisotropy at the outer radii change sharply at the time of the 
core bounce, we do not see any such sharp changes in figure 4b. These facts tell us that the 
slowing down of the increase rate of the anisotropy after the core bounce appearing in figure 
3a is mainly due to the fact that the central relaxation time becomes longer and longer as the 
cluster expands. 

In the statistical data from iV-body simulations for N = 1000 by Giersz and Spurzem (1994, 
figure 11), we can see that the anisotropy at outer Lagrangian radii reaches its maximum around 
the core bounce time, and then decreases. Such a decrease does not occur in our 2D FP models. 
Giersz and Spurzem (1994) as well as Giersz and Heggie (1994b) argued that the anisotropy in 
the outer regions is determined (at least partially) by binary activity: interactions of binaries 
with single stars, and the expulsion of stars and binaries from the core to the outer parts of 
the system. Such effects are not completely included in our models, but binaries only play the 
role of a continuous heat source. The effects of binaries on the anisotropy may be important for 
small- N systems and responsible for the fact that the anisotropy reaches its maximum around 
the core bounce in the 1000-body model. For N = 10000 clusters, there is a good agreement 
in the evolution of the anisotropy in the outer regions between the 2D FP and ./V-body models 
(Spurzem 1996; Takahashi 1996). It is not clear whether or not the anisotropy at the outer 
radii decreases after the core bounce in a 10000-body simulation (by Spurzem), because the 
simulation has not yet been continued enough beyond the core bounce. 

Figure 5 shows the evolution of the Lagrangian radii and the core radius for N = 20000. 
In this case, the post-collapse evolution of the inner Lagrangian radii is not self-similar, but 
gravothermal oscillations occur. However, the outer Lagrangian radii expand nearly self-similarly 
after the first collapse, just as in the case of N = 5000. The evolution of the anisotropy A for 
N = 20000 is shown in figure 6. The anisotropy at each of the inner Lagrangian radii reaches 
a higher value at the first collapse time for TV = 20000 than for N = 5000. This is because the 
core collapse proceeds to more advanced stages and the anisotropy penetrates into more inner 
regions (cf. figure 3 of Paper I) for N = 20000. After the first core collapse the anisotropy at 




(5) 
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the inner Lagrangian radii oscillates with the core oscillations. The anisotropy increases as the 
core contracts and decreases as the core expands. One may be interested in the fact that the 
anisotropy even at 30 and 40% radii shows a sign of oscillations, while these radii, themselves, 
show no clear sign of oscillations in figure 5. If we magnify figure 5, however, we can see that the 
radii actually oscillate with very small amplitudes. That is, we can hardly see the oscillations 
of the 30% and 40% radii in figure 5, simply because their amplitudes are too small. 

Next, we consider a gravothermal expansion phase. The mechanism of gravothermal oscil- 
lations was already clearly explained in a seminal work by Bettwieser and Sugimoto (1984). 
The key feature of the gravothermal oscillations is the appearance of a temperature inversion 
(i.e., an outward increase of the temperature) which causes a gravothermal expansion. Very 
recently, it has been clearly shown that a temperature inversion actually appears in a real N- 
body system of iV = 32000 (Makino 1996). Figures 7 and 8 show the evolution of the velocity 
dispersion (or temperature) and density profiles when a gravothermal expansion occurs in the 
case of N = 40000. At t = 18.74 t^, a temperature inversion has just appeared, and the 
gravothermal expansion has started. At t = 18.84 i r h,i the amount of the temperature inversion 
is nearly maximum. At this time, of slightly exceeds of m the region of the temperature hump. 
This is because the radial velocity dispersion at the hump reflects a lower central temperature 
more than the tangential one. The temperature inversion almost disappears at t = 19.04 t r h : i- 
This indicates the end of the gravothermal expansion, and a normal isothermal core appears 
again. 

Paper I showed that the density profile in the outer halo is approximated by a power law, 
p r -3.5 ^ c £_ Spitzer, Shapiro 1972), after the rapid development of anisotropy in the halo 
from the isotropic initial conditions. Figure 9 shows the density profiles at three epochs after 
the core collapse for N = 5000. It seems that the halo density profile further approaches the 
power law p oc r~ 3 - 5 after the collapse. As we can see in figure 2, the density profile evolves 
self-similarly well after the core collapse. (Even when gravothermal oscillations occur, the halo 
expands nearly self-similarly.) Therefore, in well-relaxed isolated clusters, the halo density profile 
is always approximated by a r -3,5 power law. Paper I also showed that the tangential velocity 
dispersion profile in the outer halo is reasonably approximated by a power law, of oc r -2 , though 
this approximation is not as good as the approximation for the density (see figures 7 and 8 of 
Paper I). This power law can be applied for post-collapse clusters as well. Actually, the velocity 
dispersion profile in the halo changes little after the collapse. 

5 Conclusions and Discussion 

In the previous paper (Paper I) an improved numerical code for solving the orbit-averaged 
FP equation in energy-angular momentum space was developed in order to study the evolu- 
tion of star clusters which have anisotropic velocity distributions. Numerical simulations were 
performed by using the code for the pre-collapse evolution of single-mass globular clusters in 
Paper I. In this paper, following Paper I, the post-collapse evolution of single-mass clusters was 
considered. The effect of three-body binaries was incorporated in the code as a heat source. 
Actually, two partially different codes were developed in Paper I. They differ in the scheme 
for solving the FP equation, itself: one uses the finite-difference scheme and the other does the 
finite-element scheme. The two codes have similar numerical accuracy for pre-collapse calcula- 
tions. For post-collapse calculations, however, the finite-element code is worse concerning energy 
conservation than the finite-difference code. Although this difficulty of the finite-element code 
may be removed by using higher-order basis functions, such efforts were not made in this study. 
By using the finite-difference code we can perform post-collapse calculations with reasonable 
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numerical accuracy (within a few percent error in energy conservation). 

There is no significant difference in the evolution of the central density between the ID and 
2D FP models as far as we studied for N = 5000, 10000, 20000, and 40000. (However, there is a 
difference in the core collapse time, as described in Paper I; the collapse time is a little longer in 
the 2D model.) In particular, the qualitative features of gravothermal oscillations are common 
to the ID and 2D models. The appearance of a temperature inversion in the 2D model is similar 
to that in the ID model. However, a slight anisotropy appears in the region of the temperature 
hump: the tangential velocity dispersion exceeds the radial one. This is an opposite anisotropy 
to the usual one, and a consequence of the lower central temperature. The opposite anisotropy 
disappears along with the disappearance of the temperature inversion. 

Clusters expand nearly self-similarly as a whole well after a core collapse. In fact, the 
expansion is consistent with the self-similar expansion law, r oc i 2 / 3 . The core-halo structure is 
more developed in the 2D model than in the ID model. However, the evolution of the half-mass 
radius in the two models is almost identical, if the time of one model is scaled so that the collapse 
times in the two models should coincide. The density profile in the outer halo is approximated 
by a power law, p oc r -3 ' 5 , after the core collapse as well as before it. 

The anisotropy at the inner Lagrangian radii decreases during a rapid core-expansion phase 
just after the core bounce. When the core expansion is stable (e.g. for N = 5000), the anisotropy 
at each of the inner radii settles to a roughly constant value, because the inner radii expand 
self-similarly as well as the half-mass and outer radii. When gravothermal oscillations occur (e.g. 
for N = 20000, 40000), the anisotropy at the inner radii oscillates with the core oscillations. 
The anisotropy at the outer Lagrangian radii continues to increase slowly after the core bounce, 
whether the core oscillations occur or not. The rate of anisotropy increase at the outer radii slows 
down as the cluster expands. This is mainly because the central relaxation time gets longer. If 
we measure time in units of the central relaxation time [see equation (H)], the anisotropy increase 
rate at the outer radii is almost constant, except for the initial epochs of the calculations when 
the anisotropy increases very rapidly. 

There are other currently-working numerical codes which can deal with the velocity anisotropy: 
one of them is Spurzem's code, which is based on the anisotropic gaseous model (e.g. Spurzem 
1996); the other is Giersz's code, which solves the FP equation by a Monte-Carlo technique 
(Giersz 1996). On the other hand, Giersz and Heggie (1994a, b) showed that the combination 
of a large number of iV-body simulations leads to results of high statistical quality which can 
give valuable information concerning the theory of stellar dynamics. Some comparisons between 
iV-body, isotropic/anisotropic gaseous, and isotropic FP models were made for isolated one- or 
two-component clusters (Giersz, Heggie 1994a, b; Giersz, Spurzem 1994; Spurzem, Takahashi 
1995). Those comparisons showed that the results of the FP and gaseous models are generally in 
good agreement with the statistical data of iV-body simulations. However, differences between 
the statistical models and the iV-body models remain in some other respects. 2D FP models 
may give a better agreement with iV-body models. Comparisons of the 2D FP models with 
the anisotropic gaseous and TV-body models are now in progress. A preliminary result of such 
comparisons for iV = 10000 models is presented by Spurzem (1996) and Takahashi (1996). For 
example, concerning the evolution of anisotropy in the halo, the agreement between the 2D FP 
and iV-body models is very good. This fact supports the reliability of our 2D FP models. 

Through Paper I and this paper, we have investigated the evolution of realistic anisotropic 
models of globular clusters. However, these models are unrealistic in some respects. They 
do not consider the stellar-mass spectrum, the galactic tidal field, and the effects of tidal and 
primordial binaries. The stellar mass-loss may also have an important effect on the initial 
evolutionary stages of clusters (Chernoff, Weinberg 1990). More realistic models incorporating 
some or all of these various effects will be studied in the future. 
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(a) N=5000 
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(b) N=10000 
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Figure 1: Evolution of the central density for (a) N = 5000, (b) N = 10000, (c) (c) N = 20000, 
and (d) N = 40000. The solid curves are the results of the 2D FP calculations, and the dotted 
curves are the results of the ID FP calculations. The time is measured in units of the initial 
half- mass relaxation time i r h,i- 
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Figure 2: (a) Evolution of Lagrangian radii containing 1, 2, 5, 10, 20, 30, 40, 50, 75, and 90% 
of the total mass, for N = 5000. The solid curves are the result of the 2D FP calculation, and 
the dotted curves are the result of the ID FP calculation. The core radii are also plotted by 
the dashed curve (2D) and the dash-dotted curve (ID), (b) Same as (a), but the time of the ID 
calculation is scaled so that the collapse time in the ID calculation should coincide with that in 
the 2D calculation. 
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Figure 3: Evolution of the anisotropy parameter, A = 2 — 2a^/a^, at the (a) inner (1-20%) 
and (b) outer (30-90%) Lagrangian radii, for N = 5000. 
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Figure 4: Same as figure 3, but the abscissa is the elapsed number of actual central relaxation 
times, r [see equation (5)]. The core bounce occurs at r « 2400. 
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Figure 5: Evolution of the Lagrangian radii in the 2D model for N = 20000. The dashed curve 
represents the core radius. 



16 



i I r 



i i i r 



.35 



b -3 



.25 



18.74 



18.84 



19.04 




J I I I 



-3 



-1 







log r 



Figure 7: Evolution of the velocity dispersion (or temperature) profile when a temperature 
inversion appears (the times are indicated in the figure in units of the initial half-mass relaxation 
time), in the 2D FP model for N = 40000. The solid and dotted curves are the radial and 
tangential velocity dispersions, respectively. 
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Figure 8: Evolution of the density profile which corresponds to the velocity dispersion profile 
shown in figure 7. 
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Figure 9: Density profiles at three epochs after the core collapse in the 2D FP model for 
N = 5000. The dotted, dashed, and solid curves represent the profiles at t/t r h,i = 17.9 (just 
after the collapse time), 28.4, and 42.8, respectively. The asymptotic line p oc r~ 3 - 5 is shown for 
a comparison. 
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